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An  algorithm  is  presented  for  the  construction  of  conformal  mappings  from  arbitrary  simply- 
connected  regions  in  the  complex  plane  onto  the  unit  disk.  The  algorithm  is  based  on  a  combi¬ 
nation  of  the  Kerzman-Stein  integral  equation  '(see  (ltflT end  the  Fast  Multipole  Method  for  the 
evaluation  of  Cauchy-type  integrals  (see  [lfrr  0,  8»  Previously  published  methods  of  this  type 
have  an  asymptotic  CPU  time  estimate  of  the  order  0(n*),  where  n  is  the  number  of  nodes  in  the 
discretization  of  the  boundary  of  the  region  being  mapped.  The  method  we  present  has  an  estimate 
of  the  order  O(n),  making  it  an  approach  of  choice  in  many  situations.  The  performance  of  the 
algorithm  is  illustrated  by  several  numerical  examples. 
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1.  Introduction 

Conformal  mappings  have  been  a  popular  analytical  tool  in  the  theory  of  partial  differential 
equations  for  more  than  a  century  (see,  for  example,  [7,  10,  13]).  As  a  mathematical  problem,  the 
construction  of  a  conformal  mapping  from  a  more  or  less  arbitrary  simply  connected  region  onto 
a  disk  is  solved  satisfactorily  by  a  whole  range  of  methods  (see  [18,  10]).  The  connection  between 
conformal  mappings  and  Szego  and  Bergman  kernels  has  also  been  well  known  for  a  long  time.  For 
practical  purposes,  constructing  a  conformal  mapping  of  a  region  C  onto  a  disk  is  equivalent  to 
evaluating  a  column  of  the  Szego  kernel  corresponding  to  H  (see,  for  example,  [7,  10,  17]). 

Unfortunately,  constructing  conformal  mappings  (or  evaluating  Szego  kernels)  numerically  is 
by  no  means  a  trivial  problem,  and  during  the  last  several  years  the  interest  in  it  has  increased 
significantly  (see  [18,  10]).  Several  algorithms  have  been  proposed,  and  most  of  them  encounter 
one  of  the  following  two  difficulties.  Most  algorithms  based  on  the  Schwartz-Christoffel  formula 
(or  its  variants)  involve  the  iterative  solution  of  a  system  of  non-linear  equations  and  suffer  from 
the  usual  affliction  of  such  methods.  Namely,  unless  a  good  starting  point  is  chosen,  the  method  is 
likely  to  fail  to  converge.  Algorithms  based  on  linear  first  kind  integral  equations  (see  [14])  avoid 
this  problem  but  lead  to  systems  of  linear  algebraic  equations  with  high  condition  numbers,  with 
the  attendant  loss  of  precision,  deterioration  of  convergence  of  iterative  algorithms,  etc. 

The  approach  we  take  is  based  on  the  result  originally  published  in  [11],  and  is  closely  related 
to  the  subsequent  work  published  in  [19,  12].  Given  a  fairly  general  region  tl  in  C1,  a  linear 
second  kind  integral  equation  is  derived  in  [11],  whose  solution  is  one  column  of  the  Szego  kernel 
associated  with  fi.  The  numerical  solution  of  linear  second  kind  integral  equations  leads  to  linear 
systems  with  asymptotically  bounded  condition  numbers  (see  Section  3  below),  and  does  not  involve 
non-linear  iterative  schemes  that  might  or  might  not  converge.  In  [19,  12],  the  analytical  apparatus 
of  [11]  is  used  to  construct  two  versions  of  the  algorithm  for  the  numerical  evaluation  of  conformal 
mappings  of  arbitrarily  shaped  simply-connected  regions  onto  a  disk.  The  version  described  in  [12] 
requires  order  0(n3)  operations  to  construct  the  mapping,  where  n  is  the  number  of  nodes  in  the 
discretization  of  the  boundary  of  the  region,  and  the  version  described  in  [19]  requires  order  0(n2) 
operations. 

In  this  paper,  we  describe  a  modification  of  the  algorithm  of  [19,  12]  with  an  asymptotic 
operation  count  proportional  to  O(n).  The  approach  is  based  on  discretizing  the  integral  equation 
of  [11]  via  the  Nystrom  method,  and  solving  the  resulting  system  of  linear  algebraic  equations  by 
means  of  the  conjugate  residual  algorithm,  which  is  quite  similar  (but  not  identical)  to  the  scheme 
used  in  [19].  Each  iteration  of  the  process  requires  that  the  matrix  of  the  system  be  applied  to  a 
vector,  which  is  normally  an  order  n2  procedure,  since  the  matrix  in  question  is  dense.  However, 
it  turns  out  that  the  the  Fast  Multipole  Method  (see  [15,  9,  3,  8])  can  be  used  to  speed  up  the 
process,  resulting  in  an  order  O(n)  algorithm  for  the  construction  of  a  conformal  mapping  of  an 
arbitrarily  shaped  region  in  Cl  onto  a  disk.  A  computer  program  implementing  the  procedure  has 
been  written,  and  we  present  numerical  examples  demonstrating  that  the  computation  times  do 
indeed  scale  linearly  with  n,  and  that  even  large-scale  problems  result  in  acceptable  computation 
times. 

2.  Mathematical  Preliminaries 

In  this  section  we  introduce  the  mathematical  apparatus  to  be  used  in  the  remainder  of  the 
paper.  D  will  denote  the  closed  unit  disk  centered  at  the  origin.  Suppose  that  Q  is  a  simply 
connected  bounded  region  in  C 1  whose  boundary  3Q  is  smooth,  and  a  is  a  point  in  the  interior  of 
fl.  Suppose  also  that  7  :  [0,  L\  — *  R*  is  a  paramaterization  of  dCl  by  its  length.  /„  will  denote  the 
conformal  mapping  fa  :  0  — ►  C1  such  that  /„(a)  =  0,  f'a{a )  is  real  and  positive,  and  fa(Cl )  =  D. 
The  existence  and  uniqueness  of  /„  we  known  as  the  Riemann  mapping  theorem  and  can  be  found, 
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for  example,  in  [4]. 

In  accordance  with  standard  practice,  we  will  denote  by  L*(30)  the  linear  space  of  all  square- 
integrable  functions  on  30,  and  by  f/*(dCl)  the  subspace  of  LJ(30)  consisting  of  all  functions  on 
30  that  are  restrictions  of  analytic  functions  on  O. 

The  Cauchy  projector  associated  with  the  region  O  is  a  mapping  H  :  (30)  — ►  X2(dCl)  defined 

by  the  formula 

H(u(w))  =  f  H(w,z)u(z)dat,  (2.1) 

Jo  n 


where 


rr,  X  1  V(*) 

=2 


is  referred  to  as  the  Cauchy  kernel,  Y(z)  is  the  counterclockwise  unit  tangent  to  30  at  z,  do  is  the 
arc  length  on  30,  and  i'{z)dot  =  dz.  The  Cauchy  projector  is  an  oblique  one,  and  the  orthogonal 
projector  S  :  L2(3 O)  — »  M2(dCl)  is  referred  to  as  the  Szego  projector.  The  Szego  projector  is  known 
to  be  a  Cauchy  type  integral  operator,  i.e.,  it  can  be  represented  by  the  formula 


S(u(tt/))  =  f  S(w,z)u(z)dot, 
Jon 


and  S  :  O  x  O  — *  C1  is  usually  referred  to  as  the  Szego  kernel. 

The  Szego  projector  is  related  to  the  Cauchy  projector  by  the  formula 

S  +  (H*  -  H)  S  =  H* ,  (2.4) 

and  therefore  each  column  of  the  Szego  kernel  satisfies  the  second  kind  integral  equation 

S(z,a)  +  f  CH(t,z)-H(z,t))S(t,a)dat  =  H(a,z),  (2.5) 

Jteo  n 

with  z  €  dfl.  The  kernel  of  equation  (2.5)  is  known  as  the  Kerzman-Stein  kernel,  and  is  defined  by 

the  formula  _ 

A(z,w)  =  H(w,z)  -  H(z,w),  for  z  ^  w  ^  ^ 

=  0,  for  z  =  w. 

A  is  smooth,  skew-Hermitian,  and  is  identically  zero  if  0  is  a  disk. 

Evaluation  of  Szego  kernels  for  regions  in  the  complex  plane  is  closely  related  to  conformal 
mappings  of  such  regions  into  the  unit  disk,  and  following  are  the  classical  formulae  implementing 
this  connection.  For  any  o  €  0\30, 


S(o,a)=  f  S(z,  a)  S(z,  a)dox. 
Jon 


For  any  z  €  0, 


For  any  z  €  30, 


Remark  2.1. 
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Obviously,  given  the  function  'b(z)  =  S(o,z),  the  conformal  mapping  fa  :  Cl  D  can  be 
easily  obtained  either  by  quadratures  using  formulae  (2.8)  and  (2.7)  or  directly  by  combining  the 
expressions  (2.8),  (2.7)  and  (2.9).  Therefore,  the  problem  of  finding  /„  has  been  reduced  to  the 
problem  of  solving  equation  (2.5). 

The  proofs  of  all  results  in  this  section  can  be  found  in  {11]  and  [12]. 


3.  Numerical  Preliminaries 

In  this  section,  the  numerical  techniques  used  to  solve  the  integral  equation  (2.5)  and  compute 
the  mapping  /„  on  the  boundary  of  fl  are  described.  The  Nystrom  algorithm  is  used  to  discretize 
the  integral  equation,  and  the  conjugate  residual  algorithm  is  used  to  solve  the  resulting  system 
of  complex  linear  algebraic  equations.  Each  iteration  of  the  conjugate  residual  algorithm  includes 
a  matrix-vector  multiplication,  which  in  existing  algorithms  has  a  cost  proportional  to  n2.  This  is 
the  only  computation  which  forces  the  operation  count  in  the  conjugate  residual  algorithm  to  be 
proportional  to  n2  rather  than  n.  The  Fast  Multipole  Method  (see  [15,  9,  3])  provides  a  technique  for 
multiplying  Hilbert  matrices  by  arbitrary  vectors  in  time  proportional  to  n,  and  we  use  this  method 
to  rapidly  apply  the  matrix  to  the  vector  in  the  conjugate  residual  algorithm  (see  Section  3.2).  Once 
the  equation  (2.5)  is  solved  numerically,  we  obtain  the  mapping  /„  by  combining  formulae  (2.7), 
(2.8)  and  (2.9). 

3.1.  The  trapezoidal  quadrature  rule  for  periodic  functions.  We  will  define  an 
n-point  quadrature  rule  q  on  the  interval  [0,  L\  as  a  finite  sequence  of  pairs  (z,,  to,},  i  =  1, . . .  ,n, 
where  z,-  €  [0,  L)  for  all  i  €  [l,n].  For  a  function  <f>  •  [0,  L]  -*  C1,  we  will  view  the  sum 

n 

n(<t>)  =  (31) 

«=i 

as  an  approximation  to  the  integral 

f  <f>[x)dx.  (3.2) 

Jo 

The  family  of  quadrature  formulae  rjn  =  {zm,  u/m),  i  =  1, . . . ,  n,  has  a  rate  of  convergence  m 
(m  >  1)  for  the  function  <t> :  [0,  L\  — *  f?1  if  there  exist  A  >  0  and  N  >  0  such  that 

\V n (4>)  -  *(*)  dx  1  <  (3.3) 

for  all  n  >  N.  The  n-point  trapezoidal  quadrature  rule  is  defined  by  the  formulae 

(3.4) 

for  i  <=  (2,r»  -  1]. 

Using  the  Euler-MacLaurin  formula,  it  is  easy  to  show  that  if  the  function  <j> :  [0,  L\  — ►  R1  is  periodic 
and  has  k  continuous  periodic  derivatives  (k  >  1)  then  the  trapezoidal  rule  for  this  function  has 
order  of  convergence  Jk  +  1 . 

3.2.  The  Nystrom  algorithm.  The  Nystrom  algorithm  associated  with  an  n-point 
quadrature  formula  tj  =  {z,,  u>,  }  replaces  the  integral  equation 


wi  =  wn=  —  aT*d 


for  i  =  1 
Wi  =  L/N 


(3.6) 


I 


with  the  system  of  linear  algebraic  equations 

n 

+  Yi  Wi  '  ■  ^(r>)  =  0(z*)> 

3=1 

with  i  =  1, ...  ,n.  We  denote  the  matrix  of  the  discretized  system  (3.6)  by  Bn>  and  view  the  solution 
•  •  • » <t>n  of  (3.6)  as  an  approximation  to  the  solution  </>  of  (3.5)  at  the  nodes  zj, . . . ,  zn.  If  (3.5)  has 
a  unique  solution,  then  for  a  wide  class  of  quadrature  formulae  r)n  and  sufficiently  large  n,  system 

(3.6)  also  has  a  unique  solution.  Furthermore,  under  fairly  broad  assumptions,  the  convergence 
rate  of  the  Nystrom  algorithm  is  the  same  as  the  convergence  rate  of  the  quadrature  formula  which 
it  is  based  on  (see  [2]). 

Applying  the  Nystrom  algorithm  based  on  the  n-point  trapezoidal  quadrature  rule  to  the 
equation  (2.5),  we  obtain  the  system  of  equations 

<t>i  +  \  A (*i.  *i)  •  =  H(a, Zi),  (3.7) 

3=1 

with  t  =  l,...,n,  where  <f>j  is  viewed  as  an  approximation  to  S(z,,  a).  The  system  of  equations 

(3.7)  can  be  written  as  the  linear  system 

Bn*  =  Hn,a,  (3.8) 


where  tf>  is  viewed  as  the  approximation  to  one  column  of  the  Szego  kernel  (S(zi,a), . . .  ,S(zn,  a)) 
and  =  (R(a,  zx), . . . ,  H(a,  zn)). 

By  the  following  theorem  (see,  for  example,  (2]),  the  condition  number  of  the  matrix  Bn  is 
asymptotically  bounded. 

Theorem  3.1.  Suppose  that  K  :  (0,  L\  x  (0,  L|  — »  Rl  is  a  c*-function  and  that  equation  (3.5)  has 
a  unique  solution.  Suppose  further  that  the  system  of  linear  equations  (3.6)  with  corresponding 
coefficient  matrix  B„  has  been  obtained  by  applying  the  Nystrom  algorithm  based  on  the  trapezoidal 
quadrature  rule  to  (3.5).  Then 

lim  k(Bn)  =  a,  (3.9) 

T%-“ *00 

where  0  <  a  <  oo  is  some  real  number  and  k(Bn)  denotes  the  condition  number  of  the  matrix  Bn. 

3.3.  Representation  of  the  boundary  of  ft.  The  discretized  equation  (3.7)  assumes  that 
the  nodes  zi,...,z„  are  equispaced  along  the  boundary  dft.  In  our  implementation,  we  assume 
that  the  boundary  of  ft  is  specified  numerically  by  nodes  wi,. ..,  wm  for  some  m  >  0.  To  generate 
the  equispaced  nodes,  we  resample  the  curve  by  passing  periodic  splines  under  tension  through  the 
nodes  wi,...,u>m,  and  generating  nodes  zi, ...  ,z„  which  are  equispaced  in  arc  length  along  the 
splines. 

3.4.  Rapid  application  of  the  matrix  B„  to  arbitrary  vectors.  In  this  section  we 
describe  the  application  of  the  Fast  Multipole  Method  (see  [9,  3])  to  the  multiplication  of  the  matrix 
B„  of  equation  (3.8)  by  the  vector  t/'  in  time  proportional  to  n. 

For  points  zi, . . . ,  z„  €  C1 ,  a  Hilbert  matrix  H  is  defined  by  the  formula 


H(i,j)  =  — - — ,  for  »'  ^  j, 

Zi  -  z, 

=  0  for  i  =  j. 


(3.10) 


A  description  of  the  Fast  Multipole  Method  (FMM)  together  with  the  proof  of  the  following  theorem 
can  be  found  in  [3] . 
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Theorem  3.2.  Let  A  denote  a  Hilbert  matrix  of  order  n,  and  let  x  be  an  arbitrary  vector.  The 
Fast  Multipole  Method  applies  the  matrix  A  to  vector  x  in  0(n  •  log{  j))  operations  where  e  is  the 
precision  of  calculations. 

Suppose  now  that  we  would  like  to  efficiently  apply  the  matrix  Bn  of  system  (3.8)  to  a  vector. 
Due  to  (2.2)  and  (2.6),  the  equations  (3.7)  can  be  rewritten  in  the  form 

* +  fWgttr* +  <3ii> 

We  denote  the  Hilbert  matrix  associated  with  the  points  zj, . . . ,  z„  by  Un,  and  the  Hilbert  matrix 
associated  with  the  points  zj, . . . ,  zn  by  Vn.  We  denote  the  vector  {xpi  •  't'(zi),  ■ . . ,  t/>n  •  V(zn))  by  x- 
It  is  easy  to  see  from  formula  (3.11)  that 


Bn  ip  =  xp  +  ctUnrp  +  c2V„x, 


L  LI 

ci  = - ==r-  and  Ci  —  —  - — 

n  2xi  n  2xt 


(3.12) 


(3.13) 


From  equation  (3.12)  it  is  clear  that  two  applications  of  a  Hilbert  matrix  to  a  vector  are  required 
to  apply  the  matrix  Bn  to  xp.  Consequently,  the  matrix  Bn  can  be  applied  to  the  vector  xp  in  order 
O(n)  operations. 

3.5.  Corrugate  residual  algorithm.  The  conjugate  residual  method  (CR)  is  an  iterative 
method  for  solving  a  system  of  linear  equations  Ax  =  f  with  i,/6C"  and  A  a  normal  n  x  n-matrix 
(see  [5]).  We  denote  the  ith  iterate  by  *,•  and  the  residual  at  the  itK  step  by  r,-  =  y  -  Ax,.  At  the  ith 
step  the  CR  algorithm  minimizes  the  residual  ||r<j|j  over  the  t-dimensional  subspace  of  Cn  spanned 
by  the  vectors  r0,  Aro, . . ,  A’_1r0. 

From  equation  (3.7),  it  is  clear  that  the  matrix  Bn  of  system  (3.8)  is  normal  since  it  is  the  sum 
of  the  identity  matrix  and  a  skew-Hermitian  matrix,  and  hence  the  conjugate  residual  algorithm 
can  be  used  to  solve  (3.8)  (see  [5]).  Moreover,  if  Eh(xi)  denotes  the  L2 -error  at  the  tth  iterate,  it  is 
known  [6]  that 

5  2IttS1<£i(io)-  (3  u) 

Combining  (3.14)  with  Theorem  3.1,  it  is  easy  to  see  that  the  convergence  rate  of  the  conjugate 
residual  algorithm  applied  to  (3.8)  is  independent  of  n. 

4.  Description  of  the  algorithm 

Following  is  a  description  of  the  algorithm  for  the  numerical  evaluation  of  the  mapping  fa  : 

an  —  dD. 

•  Resample  an  into  equispaced  nodes  using  periodic  splines  unter  tension  as  the  basis  for  inter¬ 
polation. 

•  Apply  the  conjugate  residual  algorithm  to  solve  the  linear  system  Bnxp  =  HnA  which  arises 
from  the  discretized  system  (3.7). 

•  In  the  conjugate  residual  routine,  apply  the  matrix  Bn  to  the  vector  using  formulae  (3.12)  and 
(3.13)  where  the  Fast  Multipole  Method  is  used  to  compute  the  Hilbert  matrix-vector  products 
Unxp  and  Vnx- 

•  Use  formula  (2.8)  to  compute  f'a  and  a  combination  of  formula  (2.9)  and  recursive  application 
of  the  modified  trapezoidal  rule  (see  Section  5)  to  compute  fa. 
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5.  Details  of  the  Implementation 

For  a  wide  class  of  regions  ft,  the  conformal  mapping  fa  :  3ft  — »  dD  compresses  much  of  the 
boundary  of  ft  into  very  small  segments  of  the  unit  circle.  This  phenomenon  is  known  as  crowding, 
and  is  an  inherent  property  of  the  conformed  mapping.  In  areas  of  crowding,  the  derivatives  /' 
can  be  very  small,  which  has  serious  numericed  implications.  Where  crowding  is  significant,  the 
derivative  -jr  of  the  mapping  f~l  :  dD  —*■  3ft  is  sufficiently  large  that  f~l  cannot  be  numerically 
computed  from  /„  and  f'a. 

However,  the  numerical  computation  of  /„  from  /'  is  possible.  Once  an  approximation  to  the 
Szego  kernel  has  been  computed,  the  derivative  /'  and  /„  :  3ft  — *  dD  can  generally  be  computed 
using  equations  (2.8)  and  (2.9).  However,  suppose  that  for  some  integer  t  <  n,  z,  €  3ft  is  a  point  in 
a  region  of  crowding.  Then  the  analytic  values  of  /'  are  very  small  relative  to  values  of  f'a  elsewhere 
on  the  boundary  (see  Section  6).  Therefore,  the  numerical  approximations  of  S(zi,a)  and  /'(z,) 
are  likely  to  have  a  high  relative  error,  rendering  formula  (2.9)  useless  (see  Section  7). 

Clearly,  an  alternate  method  of  computing  fa  is  integration  along  the  boundary  3ft,  and  we 
use  the  modified  trapezoidal  rule  with  4th  order  corrections  at  the  endpoints  (see  [1])  given  by  the 
formula 


f  f(x)dx  =  /(*,) 

Jx  i  •_* 


-  -  «/(*>)  +  3/(xi))  (51) 

-  -  4/(xn-l)  +  3/(x »)). 

By  separately  storing  the  right  endpoint  correction  and  the  remaining  sum  in  formula  (5.1),  the 
calculation  of  the  integral  from  z\  to  Zk+i  can  make  use  of  the  calculation  of  the  integral  from  z\ 
to  a*,  and  therefore  the  modified  trapezoidal  rule  can  be  implemented  recursively,  resulting  in  an 
order  O(n)  procedure. 

Other  than  in  areas  of  crowding,  the  formula  (2.9)  can  be  expected  to  be  more  accurate  than 
the  trapezoidal  rule.  We  combine  the  two  methods,  using  formula  (2.9)  wherever  the  derivative  is 
above  a  threshold  value,  and  recursive  integration  via  the  modified  trapezoidal  rule  otherwise. 

6.  Numerical  Results 

A  FORTRAN  program  has  been  written  implementing  the  algorithm  of  this  paper.  It  has 
been  tested  on  a  variety  of  regions,  and  in  this  section,  we  present  the  results  of  six  such  numerical 
experiments.  Tables  1-6  contain  the  quantitative  information  pertaining  to  the  Examples  1-6 
respectively.  The  first  column  of  each  table  contains  the  number  n  of  nodes  in  the  discretization 
of  the  boundary  of  the  region.  The  second  column  contains  the  CPU  time  the  algorithm  took 
on  a  VAX-8600  computer.  The  third  column  contains  the  number  of  iterations  the  conjugate 
residual  process  took  to  achieve  5-digit  precision.  The  fourth  and  fifth  columns  contain  the  relative 
accuracies  in  L2  and  L°°  norms  respectively.  Finally,  the  sixth  column  of  each  table  contains  the 
ratios 

minzgr  1  f'(z)  |  r  , 

maxj,er  I  /'(*)  I  1  } 

as  determined  by  the  calculation  with  the  corresponding  n. 


Remark  6.1.  Since  no  exact  values  are  available  for  any  of  the  conformal  mappings  con¬ 
structed  in  the  Examples  1-6  below,  evaluating  the  errors  and  convergence  rates  presents  a  certain 


difficulty.  We  have  used  a  standard  remedy,  and  estimated  the  errors  (both  L 2  and  L°°)  by  the 
difference  between  two  successive  discretizations  of  the  same  problem. 

For  each  example,  we  present  a  set  of  figures  depicting  the  equispaced  nodes  on  the  boundary 
of  the  original  region  (left  column  of  Figures  1-6)  and  their  images  on  the  circle  of  radius  1  (right 
column  of  Figures  1-6).  Following  is  a  more  detailed  description  of  each  of  the  examples. 

Example  1.  An  ellipse  with  the  aspect  ratio  2:1.  The  algorithm  has  been  applied  to  this 
region  with  n  =  16,32,. . .  ,512.  The  results  for  this  set  of  experiments  are  summarized  in  Table  1, 
and  are  illustrated  in  Figure  1. 

Example  2.  An  ellipse  with  the  aspect  ratio  3:1.  The  algorithm  has  been  applied  to  this 
region  with  n  =  16,32, . . .  ,512.  The  results  for  this  set  of  experiments  are  summarized  in  Table  2, 
and  are  illustrated  in  Figure  2. 

Example  3.  An  ellipse  with  the  aspect  ratio  10:1.  The  algorithm  has  been  applied  to  this 
region  with  n  =  16,32, . . .  ,512.  The  results  for  this  set  of  experiments  are  summarized  in  Table  3, 
and  are  illustrated  in  Figure  3. 

Example  4.  A  snake-shaped  region.  The  number  of  nodes  required  to  achieve  reasonable 
accuracy  is  larger  than  that  of  the  ellipse  examples  above.  Results  for  the  snake-shaped  region 
are  presented  for  n  =  128, 256, . . . ,  4096.  The  data  for  this  region  are  presented  in  Table  4  and  the 
corresponding  plots  are  presented  in  Figure  4. 

Exr  v  jle  5.  A  spiral-shaped  region.  As  in  Example  4,  many  nodes  along  the  bound¬ 
ary  are  required  to  achieve  reasonable  accuracy.  In  this  example,  results  are  presented  for  n  = 

200. 400. 800. . ..,  6400.  The  data  for  this  experiment  are  given  in  Table  5,  and  the  corresponding 
plots  are  presented  in  Figure  5. 

Example  6.  A  rounded  square.  The  algorithm  has  been  applied  to  this  region  with  n  = 

16. 32. .  . . ,  512.  The  data  for  this  example  are  presented  in  Table  6,  and  the  plots  are  presented  in 
Figure  6. 

The  following  observations  can  be  made  from  the  Tables  1-6  and  the  Figures  1  -  12. 

1.  The  CPU  time  requirements  of  the  algorithm  grow  linearly  with  n,  though  the  growth  is 
somewhat  erratic. 

2.  The  number  of  iterations  of  the  conjugate  residual  process  is  virtually  independent  of  the 
number  of  nodes  in  the  discretization  of  the  boundary  of  the  region. 

3.  Even  for  regions  of  fairly  complicated  shapes,  the  number  of  iterations  of  the  conjugatde 
residial  process  does  not  become  excessive. 

4.  The  algorithm  does  not  cause  substantial  round-off  errors,  and  single  precision  arithmetic 
can  be  used.  This  has  been  verified  by  repeating  the  calculations  in  double  precision.  This  obser¬ 
vation  is  typical  for  algorithms  based  on  second  kind  integral  equations  (see,  for  example,  [2],  [15], 
and  is  in  agreement  with  observations  made  in  [19,  12]. 

5.  Even  for  fairly  large-scale  calculations  (n  >  4000),  the  actual  CPU  times  required  by  the 
algorithm  are  quite  acceptable  (about  10  minutes  on  a  VAX-8600  for  n  =  4098,  and  less  than  1  hour 
for  n  =  20000).  Since  the  boundary  of  a  very  wide  range  of  regions  can  be  adequately  discretized 
by  20000  nodes,  the  algorithm  can  be  viewed  as  a  nearly  universal  tool  for  mapping  smooth  regions 
in  the  complex  plane  onto  the  unit  disk. 

7.  Crowding  and  Related  Issues 

Inspection  of  column  6  of  the  Tables  1-6  yields  a  disturbing  observation.  Namely,  even  for 
fairly  uncomplicated  regions  (such  as  an  ellipse  with  an  aspect  ratio  10:1)  the  ratio  (6.1)  tends  to 
be  extremely  small.  The  effect  is  known  in  the  literature  under  the  name  of  crowding  (see  [20]),  and 
has  received  much  less  attention  than  the  authors  feel  it  deserves.  Consider,  for  example,  Figures 
4,5.  It  is  obvious  that  nearly  all  of  the  boundary  of  the  “spiral”  is  compressed  onto  an  extremely 
small  part  of  the  circle.  For  example,  there  are  1600  nodes  in  the  discretization  of  the  boundary 
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of  the  “spiral”  in  the  left  column  of  Figure  5,  yet  only  42  images  of  these  nodes  can  be  discerned 
in  the  right  column  of  the  figure.  Similar  observations  can  be  made  about  the  “snake”  in  Figure  4 
and  about  the  ellipse  with  an  aspect  ratio  10:1  depicted  in  Figure  3.  Even  for  the  ellipse  with  an 
aspect  ratio  3:1  depicted  in  Figure  2,  the  ratio  (6.1)  is  of  an  order  IQ-2. 

One  interpretation  of  the  this  observation  is  that  while  the  mapping  f  :  Cl  —*  D  can  be 
constructed  numerically  for  a  wide  class  of  regions  Cl,  the  mapping  f~l  :  D  — *  Cl  is  a  numerically 
unmanageable  object  even  for  fairly  mild  Cl.  This  fact  has  serious  repercussions  for  attempts  to  use 
the  conformal  mappings  as  a  numerical  (as  opposed  to  analytical)  tool  for  the  solution  of  elliptic 
partial  differential  equations  (PDEs)  in  general  regions  in  two  dimensions.  Indeed,  a  standard 
prescription  for  using  conformal  mappings  to  solve  an  elliptic  PDE  (such  as  Laplace  or  Poisson 
equation)  in  a  two-dimensional  region  is  as  follows. 

a)  Find  a  conformal  mapping  f~x  :  D  — ►  fi  and  use  it  to  construct  an  induced  PDE  on  the 
disk  D. 

b)  Find  the  solution  <f>D  of  the  induced  PDE  on  D  via  an  appropriate  numerical  technique. 

c)  Construct  the  mapping  /  :  Q  — ►  D  and  evaluate  the  solution  <f>  of  the  original  PDE  on  fi 
via  the  formula 

4>(,x)  =  <t>d(f(  z)).  (7.1) 

Obviously,  if  the  mapping  /_1  can  not  be  constructed,  step  a)  in  the  above  procedure  becomes 
impossible. 

Remark  7.1 

In  fact,  it  is  known  that  for  a  region  fi,  the  ratio  (6.1)  decays  exponentially  as  a  function  of 
the  aspect  ratio  of  fi  (see,  for  example,  [20]).  It  is  also  known  that  the  mapping  /-1  viewed  as  a 
function  of  the  region  fi  is  not  even  continous  (see  [16]).  Therefore,  the  poor  numerical  behaviour 
of  f~l  is  not  surprising. 

Remark  7.2 

Obviously,  if  the  region  fi  has  an  aspect  ratio  close  to  1  (i.e.  it  does  not  deviate  from  the  disk 
too  much),  conformal  mappings  can  be  an  effective  numerical  device  for  dealing  with  elliptic  PDEs 
on  fi.  For  regions  with  high  aspect  ratios,  it  is  natural  to  attempt  to  use  model  regions  other  than 
the  disk  (such  as  rectangles  for  which  Fast  Poisson  Solvers  and  other  efficient  tools  are  available). 
While  it  is  easy  to  see  that  no  model  region  can  solve  the  problem  of  crowding  in  general,  it  is  clear 
that  certain  classes  of  regions  can  be  handled  in  this  manner. 

8.  Conclusions 

The  algorithm  of  the  present  paper  appears  to  be  the  most  efficient  of  presently  available 
numerical  tools  for  the  construction  of  conformal  mappings  from  complicated  regions  in  the  complex 
plane  onto  the  unit  disk.  The  algorithm  is  based  on  a  combination  of  the  approach  of  [11,  12,  19] 
with  the  Fast  Multipole  Method  (FMM)  described  in  [15,  9,  3],  and  has  an  asymptotic  CPU  time 
estimate  O(n),  where  n  is  the  number  of  nodes  in  the  discretization  of  the  boundary  of  the  region. 

On  the  other  hand,  our  experiments  strongly  indicate  that  the  phenomenon  known  as  the 
crowding  makes  the  use  of  conformal  mappings  as  a  numerical  tool  difficult  in  some  cases  and 
impossible  in  others.  It  is  unclear  to  the  authors  how  this  difficulty  can  be  overcome. 

The  authors  would  like  to  thank  Professor  R.R.  Coifman  and  Dr.  L.F.  Greengard  for  their 
interest  and  help. 
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Table  1:  A  2:1  Ellipse' 


nodes 

CPU  time 

iterations 

L2-er:ur 

L°°-error 

crowding 

32 

13.4 

8 

.32E+0 

.37E+0 

.60E-7  ! 

64 

15.3 

9 

.33E-2 

.39E-2 

.11E-8 

128 

21.3 

9 

.20E-3 

.20E-3 

.lOE-10 

256 

33.5 

9 

.12E-4 

.12E-4 

95E-9 

512 

57.8 

9 

.82E-6 

.85E-6 

94E-9 

1024 

121.7 

9 

.64E-6 

.66E-6 

.94E-9 

Table  3:  A  10:1  Ellipse 


nodes 

CPU  time 

iterations 

Ll-error 

L°°-error 

crowding 

128 

29.6 

17 

.32E+1 

.llE+1 

.20E-8 

256 

51.1 

17 

.70E+0 

.71E+0 

.17E-10 

512 

95.7 

17 

.lOE-1 

.lOE-1 

30E-14 

1024 

203.1 

17 

.51E-3 

51E-3 

.  15E-15 

2048 

369.7 

17 

.54E-4 

.54E-4 

98E-16 

4096 

652.1 

17 

.36E-15 

Table  4:  A  Snake-shaped  Region 


nodes 

CPU  time 

iterations 

L*-error 

L°°-error 

crowding 

200 

42.9 

wamm 

.11E+1 

.  1 1E-I- 1 

.26E-13 

400 

81.8 

.11E+0 

.11E+0 

.77E-14 

800 

127.3 

16 

.19E-2 

.19E-2 

.21E-14 

1600 

252.8 

16 

.12E-3 

12E-3 

.13E-14 

3200 

521.6 

16 

.68E-4 

.68E-4 

12E-14 

6400 

978.5 

16 

11E-15 

Table  5:  A  Spiral-shaped  Region 


nodes 

CPU  time 

iterations 

Ll-error 

L°°-error 

crowding 

16 

6.9 

6 

.30E-1 

.52E-1 

.86E-1 

32 

7.9 

6 

.83E-2 

14E-1 

47E-1 

64 

9.5 

6 

.17E-2 

43E-2 

28E-1 

128 

13.7 

6 

.35E-3 

47E-3 

20E-1 

256 

18.8 

6 

.12E-4 

45E-4 

21E-1 

512 

33.0 

6 

20E-1 

F 1 gur  e  1  a 
fl  2:1  Ellipse. 


F 1 gur  e  1 b 
R  2:1  Ellipse. 


Figure  2a 
A  3:1  Ellipse 
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A  3:1  Ellipse 


F 1 gur e  6c 
A  Rounded  Square 
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